NASA Contractor Report 178382 
ICASE REPORT NO. 87-66 


ICASE 


FRONTS PROPAGATING WITH CURVATURE DEPENDENT SPEED: 
ALGORITHMS BASED ON HAMILTON-JACOBI FORMULATIONS 


(NASA-CR-178382) FHONTS PROPAGATING NITH N88-10495 

CURVATURE DEPENDENT SPEED: ALGORITHMS BASED 

ON HAMILTON-JACOBI FORMULATIONS Final Report 

(NASA) 54 p Avail: NTIS HC A04/HF A0 1 Unclas 

CSCL 12A G3/59 0103588 


Stanley Osher 
James A* Sethlan 


Contract No. NASI -18107 
September 1987 


INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 


N/\S/\ 

National Aeronautics and 
Space Administration 

Langtoy RtteirehCcnttr 


H amp ton. Virginia 23665 


Fronts Propagating with Curvature Dependent Speed: 
Algorithms Based on Hamilton-Jacobi Formulations 


Stanley Osher 
Department of Mathematics 
University of California 
Los Angeles, California 90024 


James A. Sethian 
Department of Mathematics 
University of California 
Berkeley, California 94720 


We devise new numerical algorithms, called PSC algorithms, for following fronts propagating with 
curvature-dependent speed. The speed may be an arbitrary function of curvature, and the front can also be 
passively advected by an underlying flow. These algorithms approximate the equations of motion, which 
resemble Hamilton-Jacobi equations with parabolic right-hand-sides, by using techniques from the hyper- 
bolic conservation laws. Non-oscillatory schemes of various orders of accuracy are used to solve the equa- 
tions, providing methods that accurately capture the formation of sharp gradients and cusps in the moving 
fronts. The algorithms handle topological merging and breaking naturally, work in any number of space 
dimensions, and do not require that the moving surface be written as a function. The methods can be also 
used for more general Hamilton-Jacobi-type problems. We demonstrate our algorithms by computing the 
solution to a variety of surface motion problems. 


This work is supported in part by the Applied Mathematics Subprogram of the Office of Energy Research 
under contract DE-AC03-76SF00098, NSF under the National Science Foundation Mathematical Sciences 
Program, Sloan Foundation, NSF Grant No., DMs85-03294, ARO Grant. No. DAAG29-85-K-0190, 
DARPA Grant in the ACMP Program, ONR Grant N00014-86-K-0691, NASA Langley Grant NAG1-270. 

This work was also supported by the National Aeronautics and Space Administration under NASA Con 
tract No. NAS1-18107 while the first author was in residence at the Institute for Computer Applications 
in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 


i 



Fronts Propagating with Curvature Dependent Speed: 
Algorithms Based on Hamilton-Jacobi Formulations 


I. INTRODUCTION 


In a variety of physical phenomena, one wants to track the motion of a front whose speed depends on 
the local curvature. Two well-known examples are crystal growth [3,19,20,24,25,30,38] and flame propa- 
gation [6,18,22,23,37,39]. In this paper, we introduce, analyze, and utilize a collection of new numerical 
algorithms for studying such problems. These new algorithms approximate the equations of motion of pro- 
pagating fronts, which resemble Hamilton-Jacobi equations with viscosity terms. We demonstrate our algo- 
rithms by computing the solutions to a variety of surface motion problems. 

The background theory and numerical experimentation behind this approach have been developed in 
a series of papers, see [31,32,33,34]. In this paper, these ideas are coupled to the technology for the 
numerical approximation of hyperbolic conservation laws to produce algorithms which we call PSC 
schemes, for Propagation of Surfaces under Curvature. These new schemes allow one to follow the motion 
of an N-l dimensional surface in N space dimensions. The speed may be an arbitrary function of the curva- 
ture, and the front can also be passively advected by an underlying flow. The algorithms can be con- 
structed with any desired accuracy in space and time and do not require the front to remain a function. The 
methods are in a Eulerian framework; thus the number of computational elements is fixed at the outset. 
Topological merging and breaking is handled naturally, and the basic first order scheme is extremely sim- 
ple to program. 

As illustration of the wide applicability of such algorithms, consider the case of flame propagation, 
see [34]. A common model idealizes the burning flame as an infinitely thin boundary which separates 
regions of constant steady-state velocity, density, and temperature and propagates into the unbumt fluid at a 
speed dependent on the local curvature. The idea here is that cool convex fingers reaching out into the 
unbumt gas somehow propagate slower than do concave regions which are hot gases surrounding a small 
unbumt pocket At the same time, particles along the flame front undergo an increase in volume as they 
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bum, creating a jump in velocity across the flame front This discontinuity in the velocity field creates vor- 
ticity along the burning flame, which can be related to the local curvature, and this new vorticity field con- 
tributes to the advection of the propagating flame. Thus, there are at least two distinct ways in which the 
speed of the moving flame depends on the local curvature. 

Typically, there have been two types of numerical algorithms employed in the solution of such prob- 
lems. The first parameterizes the moving front by some variable and discretizes this parameterization into a 
set of marker points. The positions of the marker points are updated in time according to approximations to 
the equations of motion. Such techniques can be extremely accurate in the attempt to follow the motions of 
small perturbations. However, for large, complex motion, several problems soon occur. First, marker parti- 
cles come together in regions where the curvature of the propagating front builds, causing numerical insta- 
bility unless a regridding technique is employed. The regridding mechanism usually contains a error term 
which resembles diffusion and dominates the real effects of curvature under analysis. Secondly, such 
methods suffer from topological problems; when two regions "bum" together to form a single one, ad-hoc 
techniques to eliminate parts of the boundary are required to make the algorithm work. 

Other algorithms commonly employed fall under the category of "volume of fluid " techniques, 
which, rather than track the boundary of the propagating front, track the motion of the interior region. An 
example of this type of algorithm is SLIC [26]. In these algorithms, the interior is discretized, usually by 
employing a grid on the domain and assigning to each cell a "volume fraction" corresponding to the 
amount of interior fluid currently located in that cell. An advantage of such techniques is that no new com- 
putational elements are required as the calculation progresses (unlike the parameterization methods), and 
complicated topological boundaries are easily handled, see [4,32]. Unfortunately, it is difficult to calculate 
the curvature of the front from such a representation of the boundary. 

The central idea in this paper is the formulation of the correct equation of motion for a front pro- 
pagating with curvature-dependent speed. This equation in an initial-value Hamilton- Jacobi equation with 
right-hand-side that depends on curvature effects. The limit of the right-hand-side as the curvature effects 
go to zero is an eikonal equation with an associated entropy condition. By viewing the surface as a level 
set, topological complexities and changes in the moving front are handled naturally. With these equations 
as a basis, any number of numerical algorithms may be devised of arbitrary degree of accuracy, using the 
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technology developed for the solution of hyperbolic conservation laws. In particular, algorithms can be 
devised to have the correct limiting entropy-satisfying solution. In fact, some previous algorithms may be 
viewed as less sophisticated approximations to our equations of motion. 

The evolution of this approach is somewhat interesting. Motivated by the use of SLIC [26] in a 
Huyghen’s principle flame propagation scheme [4], in [31] an entropy condition was formulated for mov- 
ing fronts. In [31], it was then shown that the Huyghen’s approach was an approximation to the eikonal 
equation, which is a constant coefficient Hamilton-Jacobi equation with zero right-hand-side, and that the 
postulated entropy condition occurs naturally in this equation. Viewed from the eikonal framework, the 
inherent instability of marker particles was shown and demonstrated, see [31,34]. We then studied the 
effects of curvature on a propagating front and showed in [32,33,35] that curvature added a parabolic 
right-hand-side to the Hamilton-Jacobi equations of motion. Numerical evidence was given in [32] show- 
ing that the entropy condition formulated in [31] picked out the correct viscous limit as the curvature 
effects vanished. Attempts to approximate the solution to these equations using Lax-Friedrichs were satis- 
factory, however, the use of centered differences created spurious boundary conditions. This then led 
naturally to the higher-dimensional formulation and introduction of the higher-order upwind schemes 
employed here. 

The outline of this paper is as follows. In Section II, we give the equations of motion for propagating 
curves and surfaces in a form appropriate for numerical discretization. We then describe some past work, 
provide new proofs of some previous results, and present some new work. In Section HI, we give back- 
ground for the numerical methods for hyperbolic conservation schemes to be used and show how they can 
be used to provide solutions to Hamilton-Jacobi equations. In Sections IV and V, we use these techniques 
to approximate solutions to a variety of problems involving propagating curves and surfaces. In Appendix 
A, we discuss the inherent difficulty (linear ill-posedness) that any marker particle discretization (without 
regridding) must encounter. In Appendix B, we construct the essential non-oscillatory interpolant used in 
high order accurate approximation for general Hamilton-Jacobi equations. 


H. PRELIMINARY ANALYSIS 
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We present the equations of motion and some theoretical results about curves and surfaces moving 
with curvature dependent speed. We follow the analysis in [32] and begin with a simple, smooth, closed 
initial curve y(0) in R 2 . Let y(f ) be the one parameter family of curves, where t e [0,°°) is time, generated by 
moving the initial curve along the normal vector field with speed F , where F is a function of the curvature 
K, Let X (s ,( )=(jc (s,t),y(s ,t )) be the position vector which parameterizes yit) by s, 0 <s<S, 
X (0,t)=X ( S ,t). The curve is parameterized so that the interior is on the left in the direction of increasing j . 
With K (s ,t ) as the curvature at X (s ,/ ), the equations of motion can be written as 

* =F{K) (*/+W A y ' = ~ F(K) W+ys 1 )* i2A) 

to be solved for fe[0,°°) with X (s ,0)=^f(0), s e [0^ ] given. Here, the curvature K is defined to be 

K = ys ^ — Xs fr S Given the mapping from [0,5 ]x [0,°°) to R 2 generated by the moving curve, there 
\ x s + ys ) 

exists near t =0 an inverse mapping function / defined by t=f {x,y). The curvature K can be written in 
I terms of this function / as 


Our first result is 




K = ( 


W 


irr 


)• 


PROPOSITION 2.1 / satisfies the partial differential equation 


F 2 (f x 2 +fy 2 )= 1 

as long as the curve y stays smooth and non-intersecting. 


( 2 . 2 ) 


PROOF: The Jacobian of the mapping defined through Eqn.(2.1) is 


J = 


x t y, 
x s y. 


= F(K)(x s 2 +y 2 )* 


where K is the curvature in Cartesian coordinates. As long as this map stays smooth and one to one, we 
have f 2 + fy = t} + 1 2 = y s !J 2 + x s 2 IJ 2 = 1 IF 2 , which completes the proof. 


i 
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We notice that Eqn. (2.2) is, in general, a second order nonlinear partial differential equation to be 
solved in (x,yj) space near (x o (.s),y o (s),0), yet we are only given initial data f (x 0 (s),y 0 (s )) = constant 
on the initial curve and no information about the normal derivative of / on this curve. This seemingly para- 
doxical situation is resolved by 


PROPOSITION 2.2 Given a constant t 0 , let y 0 be the level curve of / , i.e., 

to = f (x 0 (s),y 0 (.s)). 

Then y„(s) is a characteristic curve for Eqn. (2.2). 


PROOF: Differentiating / along y with respect to s gives f x x, + f y y s = 0. Differentiating with respect to 
s again yields fxx(x s 2 ) + 2fxyX s y s + fyyfys 2 ) + f%x ss +f y y S s = 0. Using the former, we may write 
x s = a :f y y s = -ctfx for some a(^)^0, which, using the latter, gives us 

K = _ (fxXss +f )l yss) (j-2 + f 2)1/2 

Thus, the required second derivatives of / are uniquely determined on the curve from f x and f y , which in 
turn are obtained uniquely from the above. This completes the proof. 


Following [32] we define the metric g(s,t)=(x s 2 + yp) m and the angle 0 = tan -1 (y* /*,,). A simple cal- 
culation gives us 9, =gK. We differentiate Eqn. (2.1) with respect to s and rewrite the resulting system, 
using g and 0, as 

g,=QsF(Q s /g) (2.3) 

6 "Ti F( f ) ' <2 ' 4) 

Define the variation of the front at time t by 


Var(f) = |lA’(s,f)l g(s,t)ds =|l0, I ds. 
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Using this formulation, we generalize a result that first appeared in [31]. 


PROPOSITION (2.3) Consider a curve moving with speed F(K) via Eqn. (2.1). Assume F '(0)<0 and 0 
remains in the class £V'[[0,S]x[0,7’]] for 0 <f <7\ Then Var(t ) < 0. 


PROOF: The idea of the proof for smooth functions goes back to Oleinik [28] and was generalized by 
Kruz’kov [17] to the present class of functions. We shall mimic Oleinik’s proof only - the more general 
BV case follows as in [17]. Let H(s,t)= 1 if 0 s (.y,f )>0, -1 if 0 s (s ,t)<0, and 0 if 0 ,( 5 ,0=0. Then, 

U'o- 1 * = i j e -" * * 4~m F <T )H *• 

Let [$, ,s 1+ i] be an interval on which 0, >0 with 0* vanishing at the end points. Then 


• 0* 


. // 0J \ $s 


- // 0.T \ 0.V. 


F(^-)Hds = F'(— ) ^-0, I -F'(— ) 

3T 8 K 8 J 7 v 8 g" 


I . 


f 3 i a 

The first term on the right vanishes because 0, = 0 at each end point; the second term is non-positive 
because -F'(0)>0 and 0 „(j 1+] )<(K0„(5 ( ). A similar argument works on intervals for which 0,<O in the 
interior and 0, vanishes at the end points. This completes the proof. 


Using Eqns. (2.3) and (2.4), we have 

K > = - [-JjF(AT)/gL y -K 2 F(K). 
Letting F = 1 -zK , we obtain, as in [32] 


and, for e=0, 


K, = zK ss + e/s: 3 - K 2 . 


(2.5) 


fC(r t)— , 0 ) 

This becomes infinite in finite time if K (s ,0) is anywhere negative and is analogous to shock formation 
experienced in the single scalar convex conservation law. More precisely, consider the "viscous" 
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conservation law with G concave, namely 

u, +[G(u)] x =eu xx . (2.6) 

If we take e=0 (the shock case), weak solutions to 

u, + [G(u)] x =0 (2.7) 

u(x, 0) = u o (x) 

are not unique, and an additional entropy condition is needed to select the correct viscosity limit. In order 
to assure that the solution to Eqn. (2.7) be the unique limit as e— >0 of Eqn. (2.6), any of an equivalent class 
of entropy conditions is imposed [17,21,28]. The relevant one for our purposes is geometric, namely that 
characteristics flow into a shock in the direction of increasing time. This means, for a piecewise continuous 
weak solution u(x,t) having a jump moving with 

■ar - 5(0 apt — 

that G \ui) > S > G '(u r ). 

For the moving curve problem Eqn. (2.1), or equivalently, Eqn. (2.2), with F=l, we need an entropy 
condition which yields the unique limit solution as e— >0, of the problem with F=l-eK . Imagine the curve y 
as a flame separating a burnt region on the inside from an unbumt region on the outside, and an indicator 
function for the burnt region was defined to be <p(jc ,y ,t )=1 if the particle at (x ,y ) is burnt at time t , and zero 
otherwise. In [31] Sethian suggested the following entropy condition: if 4>(x,y ,f*)=l, then 4>(x,y ,< )=1 for 
t>t * , i.e., once a particle is burnt it remains burnt. This was shown to be equivalent to requiring that igni- 
tion curves flow into comers. In [32] numerical evidence was provided to show that the weak solution gen- 
erated by this entropy condition is indeed the correct limiting solution. We now prove that this is so. 

We consider a small section of the curve t=f(x,y), which, without loss of generality, we can write 
as y=Y(x ,t). We insert this into the expression fPf^=\IF 2 , arriving at 

Y, = (1+Y X 2 ) 1/2 (l+ (1 ^| )3/2 ). (2.8) 

Here, we have also chosen a positive square root Letting u = Y x and taking the x derivative of the above, 
we have [33] 

u,+[G(u)] x (2.9) 

for G(m)=-(1+u 2 ) 1/2 , G(u) concave. The criterion for the inviscid limit problem given in [31] is easily 
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seen to be that characteristics propagate into shocks for Eqn. (2.9) with e=0, that is, into corners for Eqn. 
(2.8) with F= 1. For concave / (u), this is well known to be equivalent to the statement that limits of solu- 
tions to Eqn. (2.9) (and thus Eqn. (2.8)) converge to solutions satisfying this criterion [21]. 

We may rewrite Eqn. (2.8) in the following form, namely 

Y( ~ [1+ ( i + £ V > ]([1+Y - 2]1/2) - 0 

which is a Hamilton-Jacobi equation with second-order viscosity. For the more general case of a speed 
function F ( K ), we have 


Y, " F [ (1+ vfp Kl + Y 2 ) 1 ' 2 = 0 (2.10) 

which is also a Hamilton-Jacobi equation with second-order perturbation if the speed function F satisfies 
1) F (0)*0 and 2)F'(0)*0. 

This formulation can be used to devise a numerical algorithm to approximate the solution of a curve 
propagating with curvature-dependent speed, as long as the front remains a function, using the advanced 
technology for shock dynamics. However, there is a different formulation of the problem which yields a 
different Hamilton-Jacobi type equation and does not require that the propagating front remain a function. 
Define a Lipschitz continuous function <(>(jc,jy ,r) so that at t=0, <f>(jc ,y ,0)> 1 inside the burnt region £2 , i.e., 
the region bounded by y(0), <(>(x ,y ,0)<1 outside £2, and <\>(x ,y ,0)=1 on 3 £2. Next, let 4>(x ,y ,f ) be defined by 
<j)(x ,y,t)sC where t=f (x ,y ) is defined from Eqn. (2.2) for any fixed constant C . This yields 


f - f - 

fx ~-$r fy ~~$ 7 

and hence F 2 ( + §y) = <t>, 2 . Thus, choosing the direction of propagation to be outwards, we have 


4>, -F(K) (<}>*+ 4> y 2 ) 1/2 = 0 


( 2 . 11 ) 


where 


F =F(K) = F( - 


This is also a Hamilton-Jacobi equation with second order right-hand-side. However, this different formu- 
lation allows us to compute the solution even when the front is not a function and when two "burnt" regions 
merge together. 
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Using this formulation, and the recent theory of viscosity solutions to Hamilton-Jacobi equations, 
Barles [1] has proven that the entropy condition in [31] picks out the unique viscosity solution even when 
the front is not a function. He defines <J>(x ,y ,0) = ( 1 -d (x ,y ;fl)) + + d(x,y ,Q C ) where x + =max(x ,0) and £2 C 
is the complement of Q and evolves <(> according to Eqn. (2.11) in the special case F= 1. He then chooses 
the unique viscosity solution which is characterized by the entropy condition of Crandall-Lions [5] and 
shows that the resulting surface y(t )=5Q, , defined by 3Q, = (x,y,l<t>(x,y ,0=1), evolves according to the 
entropy condition in [31]. 

Our results easily extend to initial surfaces. Suppose the surface y(0) =(x(s ijj),y(s ij 2 )>z (s 1J2)), 
moves along its normal vector field with speed F(K), yielding tOi,J 2.0=7(0- It may be rewritten as 
t=f (x,y,z), where F 2 (f 2 +f y 2 +f 2 ) = 1. Here.F =F(K), where 


1 

(ff+fS+f , 2 ?' 2 


det 


f xx fxy fxz 
fyx fyy fyi 
fix fxy fxz 


if we use the Gaussian curvature, and 


K _ (fxxify 2 HxWyy(fM?)+f»(fMy 2 y2fxyfxfy-2fxxfxfz-2fy,fyfx) 

K (fx^f?ifxW 

if we choose the mean curvature. Following the previous discussion for the propagating curve, we may 
focus on a small section of the initial surface and produce an evolution equation of the form 

Y ,-F (K) (1 + Yj? + Y y 2 ) 1/2 = 0 (2.12) 

where z=Y(x,y ,0 and K is the chosen form of the curvature. At the same time, we may once again view 

the initial surface as a level set of the function <|>(x ,yj,t) = C. More generally, to move an n -dimensional 

surface / (x ,*„)=/, we are led to the Hamilton-Jacobi-like problem 


<(>, -F(K) I V<|)l =0 

with initial data 


(2.13) 


«].(x,0) = (l-d(x,Q)) + + t/(x,f2) 

where <|>(xi ^ ,f)=l, and the curvature is chosen appropriately. Of course, it is crucial that our numeri- 

cal scheme pick out, when necessary, the correct entropy condition. 
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III. NUMERICAL METHODS 

We have seen that the problem of following a front moving with curvature-dependent speed becomes 
a Hamilton-Jacobi equation with second-order right hand side. Given an (n — 1) dimensional surface pro- 
pagating in R * , we have two formulations, namely 

1) Eqn. (2.12), which is a Hamilton-Jacobi type equation for Y in N=n - 1 space variables and applies 
when the front can be written as a function or 

2) Eqn. (2.13), which is a Hamilton-Jacobi type equation for 4> in N=n space variables, and applies 
regardless of whether the front can be written as a function. 


Thus, PSC algorithms, or Propagation of Surfaces under Curvature algorithms, rely on approximations of 

y, +//(Z>\ |/) = 0 (3.1) 

y(x, 0) = y e (x) 

with Dy = y x ,, .y^, where we have written the equations for the case F(K)= 1 for simplicity. In For- 

mulation 1), 

H(u «*) = -(l+«? + +u#y a (3.2) 

whereas in Formulation 2), 

//(mi ua,) = -(u?+.....+j^) 1/2 . (3.3) 

While Formulation 2) is more general, Formulation 1) requires one less dimension, and thus is less time- 

consuming from the point of view of numerical computations. In this section, we describe numerical 
methods that can be used to approximate the solution to Eqn. (3.1). First, we describe first order monotone 
methods for one dimension, followed by higher order methods. Then we present algorithms for first order 
monotone methods for several dimensions, followed by higher order schemes. We then show how these 
schemes can be used to solve the general case of speed function F (K). Initialization and boundary condi- 
tions are then discussed, followed by the extension of the algorithm to propagation plus passive advection. 


A. One Space Dimension 
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1) First order schemes for one space dimension 

In one space dimension, the technology for single conservation laws goes over almost directly. We 
differentiate Eqn. (3.1) with respect to the single space variable x and let u-\ \i x to produce 

u, + [H(u)] x = 0. (3.4) 

An algorithm to approximate the solution to the above is said to be in conservation form (that is, conserves 

u ) if it can be written in the form 

«7 +1 = uf - At I Ax (£7+1/2 - £7-1/2 )• (3.5) 

Here, the numerical flux function gj+in = 8 (uj-p+i Uj+q+i) must be Lipschitz and satisfy the consistency 

requirement g(u, u )=H ( u ). From here on, let y OP) be the exact (approximate) solution). 

A scheme is called monotone if the right-hand-side of Eqn. (3.9) is a non-decreasing function of all 
its arguments. It can be shown that conservative monotone schemes have no spurious overshoots nor wig- 
gles near discontinuities [16] and obey an entropy condition for limit solutions. In view of the link between 
the Hamilton-Jacobi equation and the conservation law equation in one space variable, we may easily adapt 
first order monotone schemes for shock equations to our problem (in fact, the theory goes over word-for- 
word). The simplest scheme is Lax-Friedrichs, which relies on a central difference approximation to £ , and 
preserves monotonicity through a second-order linear smoothing term. Unfortunately, this scheme is not 
upwind (to be described later), and this will turn out to be a critical requirement for boundaries. 

Thus, we begin with the canonical upwind monotone scheme, namely Godunov’s method [12]. A 
key aspect of this scheme is that the flux function g is the least viscous of all 2 point monotone fluxes [27]. 
In this scheme, g is constructed as follows. View the data as representing a piecewise constant 

function: 

Ua(x \t H ) m uf Xj- V2 <x<x J+v2 . (3.6) 

For At /Ax small enough, the initial value problem Eqn. (3.4) with u(x ,0]=u A (x‘,t n ) is a sequence of con- 

necteed Riemann problems; i.e., only adjacent constant states interact and thus may be solved "exactly" for 

one time step. This exact solution at time t n +i is then averaged over each cell to produce the numerical 

approximation uf +1 , i.e., 
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k/ +1 = - 3 j f u(x;At)dx. 

x J-v* 

Using the divergence theorem, the scheme can be put in conservation form with gj +m = H{u (xj+iaSP)). In 
other words, the numerical flux is the same as the physical flux applied to the exact solution of the Riemann 
problem Eqn. (3.4), with initial data 

U(X,0) = Uf X<Xj+\n U (x ,0) = Uj+l X>Xj+ V2 . 

We label this flux function gGor>(«/\M"+i )=gJ '+\/2 • This is clearly an "upwind" difference scheme in the 

sense that, if H’> 0, then gj +112 = f («>). likewise, if H’< 0, then gj+m = / («/+ 1 ). Another formulation of the 

above flux function [27] is simply that 

£god(w; ,«/+ 0 = Z/+ 1/2 min[x>+i/ 2 ^ (u )] 

where the minimum is taken over the interval [min(u; ,uy +1 ),max(uy,uy + i)) and Xj+v^S^j+i^j)- 

There are other useful upwind monotone schemes, see [27], which, while powerful in their own 
right, do not easily extend to several dimensions and thus are severely limited. We now present a new 
upwind monotone scheme (HJ), for the particular case when 

H(u)=f(u 2 ) 

with /'(«)< 0. Define 

gHJ ) = / ((min(n *,0)) 2 + (max(«jVi .O)) 2 ). (3.7) 

The advantage to this scheme is that it easily generalizes to several space dimensions (see below). 

It is a simple matter to put any of these schemes in terms of the Hamilton-Jacobi variable \|/. The 
numerical flux g approximates H . In the shock formulation, g must be differentiated (Eqns. 3.4, 3.5). 
However, H (and hence its numerical approximation g) appears directly in the Hamilton-Jacobi formula- 
tion (Eqn. 3.1), thus we may immediately write 

yn+i = xp«_ A , giDJV'P+'Vf). 

Here, g is any of the above numerical fluxes. (In the Hamilton-Jacobi context, it is natural to thus refer to g 
as the numerical Hamiltonian, and we shall now do so.) We note that the CFL condition is for gnj 
(At /Ax) \H'\ <1/2. 


2) Higher order schemes for one space dimension 
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Although monotone schemes have the desirable properties of conservation form, no spurious oscilla- 
tions and an entropy limit, unfortunately they are limited to first order and smear out most discontinuities. 
However, they do suggest other kinds of schemes of higher accuracy that retain these properties. 

One new class of higher order accurate algorithms was devised for conservation laws in [13,14]. 
They rely on an essentially non-oscillatory interpolant (and are thus called "ENO" schemes) and can be 
constructed to arbitrary high order. In fact, viewing them in the Hamilton-Jacobi framework results in sub- 
stantial simplification; thus we proceed directly to this setting. 

The idea is as follows. Consider the solution to y ( + = 0, with y given at /=/*. We integrate 

this in time from t=t H to t=t n+l for any fixed x and arrive at 

y(* ,t H+1 ) = y(x,r) - 1 H (y x (x,t n +s))ds. 

To approximate this procedure, let '¥ J n approximate the exact solution at time n At. We want to devise a 
function R M (x\'¥ n ) which approximates y(x,t n ) in regions of smoothness of y, up to 0( Ax) M+1 . More- 
over, this approximating function should be non-oscillatory even if y, is discontinuous, i.e., no new 
significant oscillations are introduced. We build the interpolants from the ground up as follows. For M=l, 
R\x; y ¥ n ) is defined to be the unique piecewise linear function connecting the points (x v , T'J), thus produc- 
ing precisely Godunov’s first order algorithm. For M =2, in each cell Xj<x<Xj+\, R 2 (x,'¥ n ) is the parabola 
passing through (x ; ;¥/•>, (xy+i.'Fj+i), and whichever point (Xj-i, 4'jLi) (on the left) or (xj+ 2 ,'i'J+ 2 ) (on the 
right) yields the smallest (in magnitude) second derivative, thus limiting oscillations. We store this choice 
and repeat inductively for more accuracy: that is, a cubic is obtained using the three points for the parabola, 
and an additional point, either just to the left or to the right, whichever yields a smaller magnitude third 
derivative. (The general M degree construction may be found in Appendix B.) This procedure creates a 
function which 1) interpolates 'F v n at Af+l consecutive points containing x, and x j+ i and 2) minimizes the 
magnitude of all derivatives, given the above constraint. It can be shown (see [14]) that, if y(x) is piece- 
wise C<? , has at most a finite number of isolated jumps in its derivatives and is smooth atx=x 0 , then for Ax 
sufficiently small, 

[(-^mO - (^) v * M (*;V)]*=*. = 0 (Ax) w-V+1 
where v=0,1,2,.jW . The global statement is also true, namely, 



14 


TV[ » (x ;v)] < TV[^ V (x )] + 0(Ax)* 

where TV is the total variation of a B.V. function as defined in [17]. To continue the algorithm, we then 
solve the initial value problem (Eqn. (3.1)), with 'Vfx ,0) = R M (x; 'P"), either exactly (as in the Godunov 
scheme) or approximately, using any other monotone approximation. To obtain the new 'P" +1 , we define 

% (xj ,s)ds + V? = - Jh Q¥ x (xj j )) ds 

where At is taken small enough so that only waves from adjacent cells interact. 

We simplify this method for our calculations. First, instead of solving the exact problem, we approxi- 
mate the Godunov flux by the simpler monotone flux gw The numerical time integration can be performed 
either by formally replacing higher time derivatives to arbitrary order by space derivatives, see [13], or by 
producing a non-oscillatory Runge-Kutta type algorithm [29] from the semi-discrete formulation 

^ 'F; =~ga od (£*"(*r. '*')• ' P)) • 

We note that for first order monotone approximations to a linear equation u,=-u x , the Hamilton- 
Jacobi and conservation law formulations yield the same schemes. However, differences occur for higher 
order methods. A second order approximation for conservation form gives 

= ~D-[Uj + ^jr m [D-Uj’D+Uj]] 
and 

-^- = -[£>_'P; + ^-mlP-D-'Vj.D.D+'Vj]] 

for Hamilton-Jacobi form (m is defined below). The first is only first order accurate near critical points 
because the term with m [ ] gives u z +0 (Ax), and the coefficient in the O(Ajc) is not Lipschitz continuous 
near these points. The second is uniformly second order accurate since the analogous term yields 
4' xt +0(Ax). We also have estimates 

^X\D+Uj\<0 -jj- ^ID+DAPj l<0. 



B. Several Dimensions 
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In general, the Hamilton-Jacobi equations in several dimensions cannot be recast as a simple system 
of conservation laws. However, we now design a new class of monotone upwind schemes for arbitrary H 
in several space dimensions. In our application here, namely the case H(\\f x A|r£), where/ 
is non-increasing in each of its arguments, we devise a particularly simple class of algorithms. 

1) First order schemes for several dimensions 

Crandall and Lions [5] have analyzed monotone (and hence first order) difference approximations to 
the Hamilton-Jacobi equations. As an introduction, we consider the two-space dimension discrete approxi- 
mation 

W = 'PyJt - A tgmVj-pJc-r (3.8) 

where p,qs,s ,>0, and gis thus a function of (p+< 7 +l)(r+s+l) variables. Consistency requires that g be 
Lipshitz continuous and that g(a, <a\b,b b) = H(a ,b). Crandall and Lions proved a rate of conver- 

gence result of O (Ar) 1/2 in the max norm. Their only example however is Lax-Friedrichs, which relys on a 
central difference formulation, and suffers from excessive diffusion. Unfortunately, in our solution of front 
propagation problems, the computational domain must be limited to a finite region, and thus far-field boun- 
dary conditions are required. A central difference scheme creates spurious waves at the boundary because 
it does not make use of the direction of propagating characteristics. In fact, the second author’s original 
attempt to solve the level surface Hamilton-Jacobi equation (Eqn. (2.13)) using Lax-Friedrichs suffered 
from just this problem, and this is what ultimately led to the introduction of upwind schemes. 

We begin by defining a new upwind first order generalization of Godunov’s scheme [12,27]. Let 

X#> = sgn[D*D£'I' /t ], x)# ) = sgn[DU)>'F / *] 

Define 

Hjk (« ) = Xji ) min (xftH ( u ,v )) (3.9) 

ve[mm(Dl'¥ Jk ,Dl'¥j k ), max(D> ¥;*,£* ¥;*)). 

Then this Hamilton-Jacobi-Godunov scheme has numerical Hamiltonian 

gGOD(DlV Jk ,DiVj k ,DlVj k ,DlVj k ) = xPmzx(xjt ) Hj k (u)). (3.10) 

ue [min(D£'Fy*,D*'F > *), min (Di% k ,D^ jk )) 

This is fully upwind, in that, if < 0, -2^- < 0, then the scheme looks in the proper direction, i.e., 
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g god= H(Dl'¥ji l ;Dl'¥jk). The same holds for the other three cases, and the numerical Hamiltonian g 
does not depend on the ordering of the operators. In fact, the scheme is monotone if 

r)// 

1 5 At /Ax \Hi I + At I Ay l// 2 1 . Near sonic points, i.e., points where or vanish, the Hamiltonian 
defined in Eqn. (3.9) becomes a bit complicated, and we resist reproducing the formula here. 

In our special case H(u,v) =/( u 2 , v 2 ), with / non-increasing in both variables, our one- 
dimensional HJ scheme is easily extended to two dimensions through 

gw =f(((min(Dl'¥ jk , 0)) 2 +(max(D^ ; * , 0)) 2 ) , ((min(D>^ , 0)) 2 + (maxdH'P,* , 0)) 2 )) (3.11) 
which is fully upwind and monotone, subject to the CFL restriction 1 > 2[-^- \H\ I + 1// 2 I ]. 

2) Higher order schemes for higher dimensions 

We extend our higher order methods to higher dimensions by using the spatially discrete temporally 
continuous formulation obtained from our one-dimensional ENO reconstruction procedure dimension by 
dimension. Thus, for example, a second order in space method is 

(3-12) 

DlV jk +fym[DyDl'1'j k ,D>Dl'i>j k ]m'i'jk--ty L m[DU)l'¥j k ,DV>l'¥j k ]] 

where m[x j]=xif lx l<ly I and m [x ,y ]=y if lx l> ly I. There is no loss of the desirable properties if m 
is defined as above except if xy <0 in which case it is taken to be zero. We shall use this definition in the 
next section , since it yields a "smoother" flux function. 

To obtain a fully discrete algorithm of the appropriate accuracy in time, we view Eqn. (3.12) as a 
nonlinear evolution operator of the type 

^ Jk =-LWJ,k) 

and employ certain Runge-Kutta type schemes (see [29] for a theoretical justification). For example, a 
second order essentially non-oscillatory Runge-Kutta algorithm is 


I 
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yjl +1 = 1/2 ^ + 1/2 ^r 1 - -^L['F" +1 , ./,*], 

which has a slightly reduced CFL restriction from the underlying monotone algorithm. 


C. General F(K) 

For general F ( K ), we have 

y,+F(K)H(V\ |/) = 0 (3.13) 

where K is the curvature and involves terms like y** , \\f yy , \|f*y . We have found that it is necessary to 

separate F(K) into a constant term (the convection part) and those terms dependent on K, that is, 

F(K) = F o + Fi(K), where F o is a constant (possibly zero) and Fi(0)=0. Equation (3.13) then becomes 

y« + (Fo) H (Vy) = -F \(K) H (Vy). (3.14) 

While the convection term //(Vy) on the left is approximated using one of our non-oscillatory upwind 

methods, all derivatives on the right, including Vy, are approximated by central differences. The reason for 
this may easily be seen from the following illustration: Consider a circular front of initial radius one mov- 
ing with speed F{K)=K. In Formulation 2, (Eqns. 3.1, 3.3), this is one of an infinite number of concentric 
level curves: those with small radii near the center have large curvatures. Since the term K //(Vy) on the 
right depends on multiplication and division by (y£+y^) 1/2 , which is very close to zero near the origin, the 
approximation to y, and y, must be the same within K and //(Vy), otherwise, large errors result. Thus, it 
is simplest to stick to central differences throughout the right hand side. One can also show linear stability 
of this semi-discrete approximation. We make this spatially discrete algorithm fully discrete using either 
just a forward Euler time discretization (first order accurate) or a higher-order Runge-Kutta procedure, as 
in [29]. Because of the "parabolic" right-hand-side, any such method will have a somewhat smaller CF1 
restriction than the "inviscid" approximation. 


D. Initialization and Singularities 

For propagating level surfaces, we initialize y by taking y(x,0) = 1 ± d 2 , where d is the distance 
from the point x to the initial surface, and the plus (minus) sign is chosen if x is inside (outside) the initial 
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surface. 

There will be points where Vy vanishes. On the right-hand-side of Eqn. (3.14), division by zero then 
occurs in the denominator of the curvature evaluation (N.B., this does not happen in the functional 
representation (Formulation 1)). Given the above initialization function, to a first approximation such 
points are surrounded by spherical level surfaces, and we may formally evaluate the limit Fi(£)Vy as the 
radius of the concentric level spheres goes to zero. If a mesh point falls exactly on a critical point of y, this 
limit is employed directly into the scheme at that point 


E. Far-Field Boundary Conditions 

In the case of a propagating function (Formulation 1, Eqn. (3.1-3.2)), if \y x (x ,0)=0 at x\ and X 2 
(using one dimension as an example), then we can employ symmetric boundary conditions at each end. 
However, in the level surface case (Formulation 2, Eqn.(3. 1,3.3)), Eqn. (3.1) must be initialized and solved 
for all of R N . Thus the computational domain must be truncated. With a non-zero convection term S in 
Eqn. (3.14), characteristics will head outwards far from the region of interest, and our upwind schemes are 
perfectly suited for these problems. In first order schemes, no far field numerical boundary conditions are 
needed for the convection term, since the schemes "look" in the right direction. Higher order schemes 
involve a choice of directions in order to remove spurious oscillations, thus we replace 
m [DiDl'¥ j ,DiD±'¥ j] with the second order term DID’L'Vj, at the right-hand far-field boundary, etc. 

However, the curvature term on the right-hand-side of Eqn. (3.14) must be treated with some care. If 
the convection term is relatively large, instabilities in this approximation will be swept out of the domain. 
However, if Fa= 0, the boundary plays a role. If the boundary is far from the initial surface, we may ima- 
gine that the level surface passing through each boundary point is almost a sphere. Thus, we use the exact 
solution to the collapsing sphere as the far-field boundary condition to the right-hand-side of <ty. 


F. Addition of Passive Advection 



19 


Suppose the propagating front is also passively advected by an underlying velocity field 
U = (ui, „un) in N -dimensional domain space. It can be shown thatEqn. (3.13) becomes 

t|/ f +F(K)H (Vy) + U ■ V\|/ = 0. (3.15) 

Here, of course, U may depend on x and t . For the numerical results in the next section, we used first order 

upwind differencing in each term 

^ = + urD+yi. 

However, when U depends on y, the front moves itself in a non-local manner, and more sophisticated 
methods are required. We will report on the extension of our algorithms to this important problem else- 
where [36]. 


IV. MOVING CURVES 

In this section, we demonstrate the versatility of our PSC algorithms applied to a variety of test prob- 
lems involving moving curves in a plane. We use the first and second order Hamilton-Jacobi schemes 
applied to both propagating functions and level curves. In all of these examples, the only input parameters 
are the initial curves, the time step At, the space step /i=l//V point , the order of the scheme and the speed 
function F (, K ). Everything else is handled automatically by the Hamilton-Jacobi formulation. 

A. F(K)=l-eK, Propagating function, Dependence on e 

First, we show the effect of the curvature term on the formation of singularities in the propagating 
fronts. Consider the initial curve y(x ,0)=cos8tlc , 0<x£l. Using these initial data, we compute the solution 
to the initial value problem (Eqn. 3.1, 3.2) with our second order Hamilton-Jacobi scheme and 

F(K)=\-eK, where AT= Thus, the "peaks" move slower than the "troughs". Periodic boundary 

conditions are employed in this scheme. In Figures la (e=0.0), lb (e=. 025), and lc (e=.l), we graph the 
position of the front at various times. There are Npo m f=\60 mesh points in the unit interval with time step 
Af=.001. In the case e=0 (Fig. la), comers form in the moving front, and these "shocks" propagate 
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upwards. In the case e=. 025 (Fig. lb), the front stays smooth due to the curvature term. In the case e=.l 
(Fig. lc), diffusion is so large that the peaks first start moving downwards (1-etf <0) before they flatten out 
enough to propagate upwards. 

B. Error analysis 

We wish to verify that our Hamilton-Jacobi schemes are indeed first and second order. We consider 
initial data \|r(j: ,0)=cos27tc , (Kx<l, F (K)=l-tK, e=.01 and compute the solution u k (x,t ) to Eqn. (3.1-3.2) 
for h =.01/2 ? , Af =.005/2? , p =0, 1 , 2, 3, 4, 5, 6. If we assume that the exact solution v(x,t) can be written as 

v(xj) = u k (x,t) + Ch* +0(h*+ l K 
then we can estimate the order R of the method by 

2 a - v ~“* 

V - U h , 2 ’ 

Since we do not know the exact solution, we approximate v by the most accurate calculation. In Table 1, 
we give the value of R found by comparing u* and u k a for various values of h at time f=1.0, using 
u* =i /6400 (p= 6) as the exact solution. Here, errors are measured using the discrete L 2 norm. The value of R 
approaches 1 for the first order scheme and approaches 2 for the second order scheme. Although the results 
appear a little better than that, we hesitate to draw any conclusions. 


Actual Order of Hamilton-J 

Facobi Schemes 


"First" Order 

"Second" Order 

h=.01/10 


R=2.2144 

h=.01/20 


R=1.8581 

h=.01/40 


R= 1.8994 

h=.01/80 

■ 

R=2.0369 

h=.01/160 

R=1.5850 

R=2.3084 
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C. F(K)=K, Propagating functions 

We consider a periodic curve y(x ,0)=sin2iu as initial data, F (K)=K, and solve Eqn. (3. 1-3.2) using 
our second-order scheme. This corresponds to a curve moving under its curvature. This problem has been 
studied extensively in [7,8,9,10] and reduces to a relatively straightforward parabolic equation. In Figure 2, 
we draw the front at various times, demonstrating that the periodic front relaxes to a straight line with 
increasing time. 

D. Level Curve, Burning out, Development of Corners 

We consider a seven-pointed star 

y(s) = (.l+(.065) sin(7-2jw)(cos(2jw),sin(2jw)) 

s €[0,1] 

as the initial curve and solve Eqns. (3,1), (3.3) with F(K)= 1, using the initialization given in Section 
III.D. The computational domain is a square centered at the origin of side length 1/2. We use 300 mesh 
points per side and a time step A/ =.0005. At any time n At , the front is plotted by passing the discrete grid 
function to a standard contour plotter and asking for the contour ^=1. The initial curve corresponds to 
the boundary of the shaded region, and the position of the front at various times is shown in Fig. 3. The 
smooth initial curve develops sharp comers which then open up as the front bums, asymptotically 
approaching a circle. 

E. Level Curve, Motion Under Curvature 

With the same initial curve as Example IV.D above, we let F(K)=K , corresponding to a front mov- 
ing in with speed equal to its curvature. It has recently been shown [10] that any non-intersecting curve 
must collapse smoothly to a circle under this motion. With /Vpoint=300, and Af=.0005, in Fig. 4a, we show 
the front at time t=.0,.01,.02,.03,.04,.05. Here, we have scaled time by a factor of 100, because the real 
front moves so quickly. In Figures 4a-4d we show the continued evolution of the surface from t=0.0 to t=.2. 
The plots show the relaxation of the peaks and troughs and the smoothing into a circle. In Figure 5, we 
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show the results of the same motion applied to a different initial curve, namely the wound spiral traced out 
by 

7 (s) = (.le<- 10 > 0 » - (. 1 — jc ( i ))/ 20 )(cos(a (s)), sin(a ( 5 ))) 
where a(s) = 25tan -1 (10y ( s )) and 

x(s) = (.l)cos(2jw>hl y (s ) = (.05)sin(2jrs )+. 1 se[0, 1]. 

With Wpointp200 and A/ =.0001, Figure 5a shows the unwrapping of the spiral from r=0 to t=0.65. In Fig- 
ures 5a-d we show the collapse to a circle and eventual disappearance at r=.295 (The surface vanishes 
when < 1 for all ij.) 

F. Level Curve, F (AT)=l-eAf , Merging and Breaking 

Using the wound spiral initial curve in the above example, Figure 6 shows the results with 
F(K)=1-eK , lVpoint=200 and A/=.0001. Figure 6a shows the initial curve as the boundary of the shaded 
region. In Fig. 6b, the spiral expands and pinches off due to the strong convection component, separating 
into two curves, one propagating outwards and one shrinking in. In Fig. 6c, the front at r=.04 is the boun- 
dary of the shaded region. The outer front expands and the inner front collapses and disappears. In Fig. 6d, 
all that remains is the outer front which asymptotically approaches a circle. 

G. Level Curve, Passive Advection and Propagation 

Finally, we solve the passive advection plus propagation equation (Eqn. 3.15) with the initial seven- 
pointed star in Example D, F(K)= 1, and 

U = (« i(* ,y ,t ). u 2 (x ,y ,f)) = ( ~y , x ) (100(x 2 + y 2 )). 

This corresponds to solid body counterclockwise rotation around the origin with tangential velocity 1 along 
the circle with radius . 1 centered at the the origin. In Figure 7, we show the expanding and spinning star at 
various times. 


V. MOVING SURFACES 
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In this section, we use PSC algorithms to compute the evolution of several two-dimensional surfaces 
in three space dimensions. We use the initializing function given in Section III.D. 

A. Propagating Function Surfaces, F(K )=1 

We evolve the initial surface 

\|/(;t,y,0) = -.25[cos(2;w)-l][cos(2jcy)-l] + 1 

according to Eqns. (3.1), (3.2), with F(K )=1, A/ =.01, N point=50 (in each direction) and periodic boundary 
conditions. This surface is flat in the boundary of the unit square centered at the origin and has a global 
minimum at (0,0). In Figure 8a, we plot the surface at various times, showing the focussing of the 
minimum into a deep dent which then opens up. The surface moves upward with unit speed, asymptotically 
approaching a flat sheet. Next, we add curvature effects to the speed function and let F(K)=\-eK, e=.l. 
The time step is reduced to At =.0001 because of stability requirements from the addition of a parabolic 
term. In Figure 8b, we plot the surface at various times. Here, the dent is greatly smoothed due to the cur- 
vature effects, and the surface becomes flat much faster, similar to the one-dimensional case (Fig. 1). 

We then consider a saddle surface, described by 

y(jt ,y ,0) = cos(2tlc ) - cos(2ny ). 

Again, we start with F(K)=1, A/=.01, Afpoim=50 and periodic boundary conditions. In Figure 9a, we plot 
the surface at various times. Here, the rising surface develops a discontinuity passing through the saddle 
point in the y coordinate direction, corresponding to a shock in the tangent vector. Adding curvature (Fig. 
9b, F(K)=l-eK , e=.l), the shock is smoothed out, and the surface smoothly approaches a flat sheet. 

Finally, we move the saddle surface purely under its own curvature ( F(K)=K , where K is the mean 
curvature). In Figure 10, we show the front at various times and show the rapid motion toward the steady 
state flat sheet 

B. Level Surface, Sphere, F(K) = 1, Mean Curvature 

We evolve the initial surface described by the sphere of radius .5. We initialize yfjjc using the dis- 
tance function, as described in III.D; in this case, the distance function from the initial surface may be 
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analytically expressed. The computational domain is a rectangular parallelpiped with lower left comer 
(-1-1-1) and upper right comer (l.,l.,l). We evolve the surface according to Eqns. (3.1), (3.3), with 
F(K)- 1, A/=.01 and Npoint=50 points per side of the domain. The solution to this problem is just a sphere 
whose radius r is increasing at the constant rate drldt- 1. In Figure 11 (Figs. 11a and lib), we show the 
front at various times. (In the display of all of the following three-dimensional plots, the axes shown are for 
orientation purposes only, and are not necessarily located at the real origin of the figure.) As expected, the 
sphere expands smoothly. Of particular interest are the three surfaces in Fig. 1 lb. As the sphere expands, 
part of its boundary intersects the edge of the computational domain. This is reflected in the slicing of the 
level surface y=l by the sides of the box. This demonstrates the advantage of an upwind formulation: since 
information flows into the boundary, the surface does not know about the boundary, and the interior points 
of the level surface proceed unharmed. 

C. Level Surface, Torus, F (K) = 1-eK 

We evolve the toroidal initial surface, described by the set of all points (x ,y z ) satisfying 

z 2 = (R o^dxW^Ri) 2 

where R $=. 5 and /?i=.05. This is a torus with main radius .5 and smaller radius .05. Again, to save labor, 
the initialization may be analytically expressed. The computational domain is a rectangular parallelpiped 
with lower left comer (-1,-1, -.8) and upper right comer (1.,1.,.8). We evolve the surface with 
F(K) = l-e£,e=.001, A/ =.01, and #^,=90 points per* andy side of the domain and the correct number 
in the z direction so that the mesh is uniform. Physically, we might think of this problem as the boundary 
of a torus separating products on the inside from reactants outside, with the burning interface propagating 
with speed K=\-tK. Here, K is the mean curvature. In Figure 12 (Figs. 12a and 12b), we plot the surface 
at various times. First, the toms bums smoothly (and reversibly) until the main radius collapses to zero. At 
that time (T=0.3), the genus goes from 1 to 0, characteristics collide, and the entropy condition is automati- 
cally invoked. The surface then looks like a sphere with deep inward spikes at the top and bottom. These 
spikes open up as the surface moves, and the surface approaches the asymptotic spheroidal shape. Again, 
when the expanding toms hits the boundaries of the computational domain, the level surface y=l is clipped 
by the edges of the box. In the final frame (T=0.8), the edge of the box slices off the top of the front. 
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revealing the smoothed inward spike. 

D. Level Surface, Sphere and Torus, F (K) = K, Mean Curvature 

Finally, we show the flow of two surfaces under their mean curvature. This problem has been studied 
extensively, see [2,1 1]. First, we evolve the initial sphere of radius 1.0 moving with speed F(K)=K. Again, 
time is scaled by a factor of 100 because the real solution goes so fast. By symmetry arguments, the evolv- 
ing surface should be a sphere of decreasing radius which eventually disappears. In Figure 13 (Figs. 13a 
and 13b), we show the collapsing sphere at various times. As easily seen, the radius decreases slowly at 
first and very quickly at the very end. The final shape shown (T=6.6) is the smallest surface that can be 
resolved on the given mesh size. 

Finally, we evolve a toroidal initial surface under its mean curvature. The inner radius is .25 and the 
outer radius is .5. We imbed the problem in a unit cube of side length 2., and use a fairly coarse mesh of 
Npomt=45 per side. In Figure 14 (Figs. 14a and 14b), we show the surface of the front at various times. For 
our particular initial surface, the torus deflates smoothly and collapses to the ring shown at T=4.1 before it 
vanishes. 


SUMMARY 

We have presented a class of algorithms, called PSC schemes, for moving surfaces under their curva- 
ture. These algorithms rely on numerically solving Hamilton-Jacobi equations with viscous terms, using 
approximation techniques from hyperbolic conservation laws. To demonstrate our techniques, we compute 
the solution to a variety of surface motion problems. We hope that this tool can be applied in several areas, 
such as flame stretch, vortex sheet rollup, Hele-Shaw cells, and crystal growth. 
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Copies of the computer program are available from the second author. All calculations were performed at 
the University of California , Berkeley and at the Lawrence Berkeley Laboratory. 
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APPENDIX A: INSTABILITY OF MARKER PARTICLES 


Here, we analyze in some detail the difficulties inherent in a marker particle discretization of any 
Hamilton-Jacobi equation, and relate this to the motion of a front moving with constant speed. 

Consider the general Hamilton-Jacobi equation 

V«=//(V*H> (A' 1 ) 

with smooth initial data vy(jcoO ),0) + y 0 (s ). The method of characteristics tells us that xg* is constant 

along curves in x space defined by 

• x( s fi)=xo(s) (A-2) 

^- = V,//'(Vx)-«(Vx) 


Following our notation in Section II, we define g = (x s 2 + y s 2 ) 112 and 0 = tan -1 where now 


x =x (s ,t ) and ,f ). This leads to a generalization of Eqns. (2.3-2.4) with F(K)=1, namely 

0 < = 0 

g, = -j\cos0)- 3 //"(tan0)j /0, 

For a curve moving with constant speed, H"( tan(0)) = -(cos0) -3 , and the system becomes linear 


(A-3) 


0 


0 0 


0 

8 

t 

1 0 

. 


8 


(A-4) 

This is a slightly ill-posed non-strictly hyperbolic system. However, the exact solution is easily seen to be 

0 = 0o(s) g= go(s ) + f0oO) (A-5) 


which "loses” a derivative. 

This linear ill-posedness is manifested by g (s ,t ) becoming zero, which corresponds to the intersec- 
tion of characteristics in the original problem (Eqn. (A-l)) at time t =r cnemin(-go(^ )/(0oC? ))"), as in Eqn. 
(2.5). 

For general concave (//"< 0) Hamilton-Jacobi equation, this occurs at 
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crit max[-//"((Yo)x) (Vo)„] 

and is hence determined by the second derivative of the initial curve, when written as y=\y(x,0). This is a 
subtle point, which we illustrate by example. 

Consider the flat initial data 

y= Oh\|/o(jO 

The true solution to a unit speed moving front with this initial data is just x=s , y=s+t. We now consider a 
sequence of initial data which converge to the above flat data as the initial space step is refined. Take the 
initial data 

yo = -As (sin(As) 1/2 ) = yo(As) = -yo(2As) = -yo(3As) 
and yoO As )=yo(0 -4) As) for all j , and 

xo(J As)=jAs cos((A s ) 1/2 ) 

defined on the grid j Ar=0,±l,±2, • • • . As As— >0, the discrete initial data converges to the flat line \|/=0 in 
the following fashion: 

y o(j As ) = 0 + O (( As ) 3/2 ) 
yoO’ As) = 0 + 0 ((As) 1/2 ) 
xrfjAs) = x +0(As) 
xoOAj) = 1 + 0(As) 

However, a marker particle numerical scheme (without regridding) blows up after 0(l/(As) 1/2 ) time steps. 
With any relationship between A t and As of the form At=(Asy,p> 1/2, the actual numerical blow up time 
is 0(A/)0(l/(As) 1/2 ) = 0(Asy~ V2 , which goes to zero as As vanishes, even though the real solution is 
finite (in fact, zero in all derivatives) for all time. 
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APPENDIX B: GENERAL ENO CONSTRUCTION 

We build R M by induction on M . In each cell Xj<x<xj+i,R l is defined as follows: 

/? i (jc;'F") = '¥f + (x-x j )D+'¥* 

= ¥/ + <x-x y )¥« 1] 

where v[xj-\i, ^y+v] denotes the usual coefficient in the Newton interpolating polynomial. We also 

define = j,k^i x = j+l. Suppose we have defined R M ~ 1 (x ; V F' 1 ) for Xj<x<Xj+\, and also k$&r l) , ktfk’ 1 '* 
(the leftmost and rightmost indices, respectively). Then we compute 

a M =V K [x k tg- l) jc k tg- 1)H ] 
b M ='V K [x k ^. 

and proceed inductively. 

If I a M I > I b M I , then in this interval 

* w-i) 

RM(x;'¥*)=R M - l (x\'¥*) + b M n (x-x k ) 

with kiL = kj&rV-l, and k%» x = A$k _l) . 

If I a M I < I b M I , then in this interval 

R M (x;'V , ) = R M ~ i (x;'V‘) + a M Tl (x-jc*) 

* = *mm 1 

with^JL =k^T 1 \ and k^ n =ktfir l) +l- 

To summarize, in each cell Xj<x<Xj+ 1, we have constructed an essentially non-oscillatory polyno- 
mial of degree M . This polynomial is the restriction to the cell of a polynomial interpolating 0P£) at M +1 
consecutive points x v , including xj and xj+u moreover, these points are chosen so that all derivatives of the 
polynomial are as small as possible in absolute value. 
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Fig. la. F(K>-1-«K, a-0.0 


Fig. lb. F(K)-1-«K, •-0.025 


0.000 0.200 0.400 0.600 0.800 1.000 0.000 0.200 0.400 0.600 0.800 1.000 



0.000 0.200 0.400 0.600 0.800 1.000 


1.500 

1.300 

1.100 

0.90C 

0.70C 

0.50C 

0.30C 

0 . 10 * 

-0.10' 

-0.30 

-0.50 


Fig. lc. F(K)-1-«K, «-0.1 All PLOTS AT T-0.0, 0.9 (0.05) 
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fig. 4*. F(K>-K T- 0.0, 0.05 (.005) Fig. 4b. F(K)«1 T-0.05, .10 (.005) 



rig. 4c. F(K)-K T-.10, .15 (.005) 


fig. 4d. T(K)-K T-.15, .20 (.005) 
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rig. 5*. r(K)-K T- 0.0, 0.065 (.01) rig. 5b. r(K)-l T-0.065, .130 (.01) 



rig. 5c. r(K)-K T-.130, .195 (.01) 


rig. Sd. r(K)-K T-.195, .295 (.01) 
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pmgikal W® 6 , ® 
OF PO° 


rig. 6a. F(K)»l-aK T-0.0 (INITIAL CURVE) Tig. 6b. f(K)-l-«K T-0.00,0.03 (.01) 



rig. 6c. r(K)-l-*K T-0.04,0.12 (.01) 


rig. 6d. F(K)-l-aK T-0.13,0.22 (.01) 
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7 a> F(K)»x ♦ Rotation T-0 . 0, . 03 ( . 01) Fig. 7b. F(K)“1 ♦ Rotation T-0.03,0.06 (.01) 



Fig. 7c. F(K)-1 ♦ Rotation T-0.06,0.09 (.01) 


Fig. 7d. F(K)-1 + Rotation T-0.09,0.12 (.01) 



41 


ORIGINAL PAGE IS 
OF POOR QUALITY 



Figure 8s: F(K)— 1 Surface at T— 0.0, 0.3, 0.6 


Figure 8b: F(K)—1- < K , < —. 1 , Surface at T«=0 0,0.3, 0.6 
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Figure 9a: F(K)— 1 SurfAc* At T— 0.0, 0.3, 0.6 Figure 9b: F(K)— 1- < K , t — .1, SurfAce At T— 0.0, 0.3, 0.6 
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Figure 12a: F(K)“1- « K , t ■» .01, Surface at T—0.0, 0.1, 0.2, 0.3 


FIGURE 12: BURNING TORUS 
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